The Approximate Solution of Linear Differential Equations 

By MARION C. GRAY and S. A. SCHELKUNOFF 

Linear differential e(|uations with variable coefficients occur in many fields of 
applied mathematics: in the theories of acoustics, elastic waves, electromagnetic 
waves in stratified media, nonuniform transmission lines, wave guides, antennas, 
wave mechanics. The "Wave Perturbation" method described in greater detail 
elsewhere' is particularly useful in those ranges of the independent variable in 
which the "WKB Approximation" is not sufficiently accurate. The present 
paper endeavors to illustrate the remarkable accuracy of this method, particu- 
larly when compared with Picard's method. 

I. Introduction 

IN A recent paper^ the approximate solution of linear differential equations 
by a wave perturbation method was described. When the method was 
applied to equations whose exact solutions were known we were greatly 
impressed by the rapidity of convergence of the successive approximations. 
Hence the purjiose of this note is to present some illustrations in the hope 
that others may be interested and may find the proposed method an im- 
provement on those now in use. 

In essence the wave perturbation metliod dates back to Liouville-, but 
in his memoires he was interested in a problem of heat conduction involving 
a non-homogeneous differential equation with homogeneous boundary' 
conditions, whereas we consider a homogeneous equation 

y" = Hx)y (1) 

with non-homogeneous initial conditions 

y(a)= l,y'{a) = (2a) 

or 

>'(a) = y'{a) = 1, (2b) 

the solution being desired in an interval a ^ x S b. Since the solution 
for any assigned initial or boundary conditions can be expressed as a linear 
combination of the solutions satisfying (2a) and (2b) we have not imposed 
any real limitation. 

II. Theory 

Comparison of the wave perturbation method with Picard's method 
(which is essentially a hnear perturbation method) is particularly instruc- 
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tive. It will be recalled that in Picard's formulation the differential equa- 
tion (1) is replaced by an integral equation 

yix) = y{a) + {x - a)y'{a) + /" F{u)y{ii){x - u) du (3) 

where y{a) and y'{a) are assigned initial values*. Writing 

U{x) = y{a) + {x ^ a)y'{a), (4) 

LS-y) = f F{u)L„.,{u)ix - u) du, n = 1, 2, 3, • • • , 

the series 

3-0^) = Ux) + L,(x) 4- L,{x) + ■ ■ . (5) 

is shown to converge to a solution of the original equation. In practical 

applications, unfortunately, it is usually found that the successive approxi- 
mations converge rather slowly unless the interval (a, b) is small. 

In the wave perturbation method we first rewrite equation (1) in the 
form 

y" = -iS'y + [^' + n^-)]y = -P'y+f(x)y, (6) 

and instead of the integral equation (3) we use 



y(x) = y(a) cos 0{x " ^) + 3 y'(^) sin &ix - a) 
+ - / f{u)yiu) sin fi{x — u) du. 

P ■'o 



(7) 



The parameter ^ is arbitrary and might be defined in various ways. We 
have found it convenient to use the definition 

so that if F{x) is negative /3 is real and our first approximation 

Waix) - y{a) cos i3{x - a) + - y'(a) sin /3(.^■ - a) (9) 

is sinusoidal. If /'' is positive ^ is imaginary and we start with an exponen- 
tial approximation. If F changes sign in (a, b) the best procedure is to 

* This is not quite the usual form of the integral equation but it is substantially 
equivalent. 
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subdivide the interval and obtain separate approximations, though this is 
not necessary if F is predominantly of one sign throughout (a, b). To (9) 
we now add the sequence 

WM ^\ I /(«)#„_,(«) sin /3(a- - u)du, (10) 

and the series ■ .^ _ 

y{x) = W,{x) + Wi{x) + W^ix) + ■ • ■ ' ] (U) 

is the desired solution. 

The flexibihty of the wave perturbation method as compared with Picard's 
linear method hes essentially in the introduction of the variable parameter |S. 
Since we make /3 depend on the length of the interval (a, b) m which a solu- 
tion is desired the approximations may be extended over much longer 
intervals than is feasible in Picard's method. If F(x) is a slowly varying 
function throughout (a, b), so that /(a;) is small, it will be found that the 
first approximation Woix) is good, and the second, Wa + H^i is generally 
adequate. 

Another choice for /3 is 

- /3 = -^ I V^^fU dx. (12) 

— a Ja 

However, the integration m (8) will often be simpler than in (12). 

Picard's method is a special case of the wave perturbation method, with 
|3 = 0. In fact, if F{x) changes sign in (o, 6), then in some cases ^ as defined 
by (8) will reduce to zero. 

If F{x) is a rapidly varying function, or if the solution is desired over an 
mfinite interval, it is usually advantageous to transform equation (1) by 
first introducing a new independent variable 

6 = I V^Hx) d^> (13) 

and then removing the first order term in the new equation by an appro- 
priate transformation of the dependent variable. 

III. Examples 

For our illustrations we have used mainly the simple equation 

y" = -xy (14) 

whose exact solution can be expressed in terms of Bessel functions of order 
± 1/3. Since the Bessel functions are oscillatory m nature it might be 
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suggested that comparison with Picard's method is weighted in our favor. 
This does not seem to be the case, as will be illustrated in example 4 where 
the exact solution is a monotonically increasing function. It has also been 
suggested that Picard's solution might be improved by starting from a better 
initial approximation, say Wo , rather than from the linear approximation 
Lo , but we have not found any marked improvement in the succeeding 
approximations (see examples 1 and 2). The various points of interest 
will be brought out in our examples, with the accompanying figures, which 
we shall now briefly describe. In each figure the heavy curve is the ac- 
curate solution while the approximations are indicated by self-explanatory 
letters. 
Example 1, Fig. I 

y" - -.17,0 ^ X ^ 2 

y(0)^ l,y(0) ^0 

Exact solution: y(x) = r(f)3"'".i:"V-i (ir'") 

(a) Wave perturbation 

1^0 = cos .T 

Wi = ~lx cos X -\- \(l + 2x - x') si 

(b) Linear perturbation 

Z-o= 1 

(c) Linear perturbation using initial sinusoidal approximation " '^ 

Lo = cos .V = ll'o 
Li = X + .1' cos .V — 2 sin x 
Example 2, Figs. 2, 3 and 4 

y" = -xy, 2 ^ X -^6 
>'(2) ^ 0, y'{2) = 1 
Exact solution: 

y{x) = ~.m2Zx"'j-i„{lx"') - .m29\x"^Ji„{lx"') 
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(a) Wave perturbation, Fig. 2 

P^o = ism2(i; - 2) 

Wi = ——- sin 2{x - 2) + — cos 2(x - 2). 

61 16 

Figure 2 exhibits rapid pulling of the successive approximate waves to- 
ward the exact even though the interval has been chosen deliberately 
unfavorable to the straight wave perturbation method [see example (c) 
and Fig. 4 for the improved treatment]. 

(b) Linear perturbation, Fig. 3 

Lq= X - 2 

L, = ^ (16 - 16x + 4x^ - X*) 

^^ ~ "63 45 "'9~9~90504" 
Using TFo instead of Lo 

Li = i - ^ + -^ sin 2(.v - 2) + i cos 2{x - 2). 

(c) Preliminary transformation of variables, Fig. 4 

Introduce 6 = fx''^ y = d~^'% 
and the modified equation is 

Then, using for simplicity ,8=1 

vo = 2~"V sin {8 - do), do = W^p 
or 

W, = (Ix)-'" sin Ux" - 2'") 

It will be seen that Wo is a very good approximation throughout the 
range (2, 6). Adding Wi obtained from 



sei 



al/6 

vi = ^^ [cos (5 + do) {Si 28 - Si28o) - sin(^ + eo){Ci2e - Ci28o)\ 
36\/2 



the accurate curve y is reproduced. 

Li Fig. 2 the third approximation could not be distinguished from the 
accurate curve though numerically the values are not identical. For 
purposes of comparison the table of numerical values (Table A) may be 
found interesting. 
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X 


y 


Wo 


1 

sir 

D 


3 

Sir 




I 

SIC 



2.0 


0. 


0. 


0. 


0. 


0. 


2.2 


.19721 


.19471 


. 19720 


.19721 


.19721 


2.4 


.37694 


.35868 


.37668 


.37694 


.37694 


2.6 


.52056 


.46602 


.51885 


.52053 


.52056 


2.8 


.61035 


.49979 


.60442 


.61020 


.61034 


3.0 


.63236 


.45465 


.61792 


.63180 


.63232 


3.2 


.57922 


.33773 


.55169 


.57781 


.57918 


3.4 


.45287 


. 16749 


.40907 


.44995 


.45726 


3.6 


.26584 


- .02919 


.20603 


.26085 


.26561 


3.8 


.04126 


-,22126 


- .02974 


.03408 


.04087 


4.0 


-.18921 


- .37840 


-.26229 


-.19800 


-.18974 


4.2 


- .38951 


- .47580 


- .45326 


-.39852 


- .39011 


4.4 


-.52506 


- .49808 


-.56889 


-.53240 


-.52557 


4.6 


-.56943 


-.44173 


-.58697 


-.57328 


-.56964 


4.8 


-.51062 


- .31563 


-.50217 


-.51000 


-.51044 


5.0 


- .35548 


-.13971 


- .32847 


- .35076 


-.35494 


5;2 


- . 1.3068 


.05827 


- .09772 


-.12364 


- . 13006 


5.4 


- . 12052 


.24706 


. 14547 


.12725 


.12080 


5.6 


.34582 


.39683 


.35200 


.34988 


.34536 


5.8 


.49485 


.48396 


.47807 


.49525 


.49347 


6.0 


.53114 


.49467 


.49467 


.52903 


.52903 



Example 3, Fig. 5 



V" = —y+ -y, 
X- 

y(l) = hy'd) = 



1 < .V< °o 



Exact solution: y{x) = sin (.-i: — 1) + - cos {x — 1) 

(a) Wave perturbation, with the initial conditions satisfied exactly 

IFo = cos (.V - 1) 

II', = 2 sin (.V - 1) - 2 cos {x + 1) {Ci 2x ~ Ci 2) 
- 2 sin (x+ 1) {Si2x - Si 2) 

(b) Wave perturbation, matching the exact solution at infinity 

H/o - sin (x - 1) 

Wi - 2 sin (x + 1) Ci 2x - 2 cos (x + \)(Si 2x - 7r/2) 

This is an example of a solution in an infinite interval, where the per- 
turbation term is not small throughout. It is interesting to note that the 
second form gives good agreement with the accurate solution in most of 
the range of integration. 
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Example 4, Fig. 6 

y" = -\- xy, ^ a: g 2 

7(0) = 1, y'{0) = 
Exact solution: 

(a) Wave perturbation 

Wo = cosh X 

Wi= -J cosh x-\-i(x^ - 2x+ 1) sinh x 



(b) Linear perturbation 



Ld= 1 
^' -6- 



This is an example in which the exact solution is non-oscillatory yet even 
in the short interval (0, 2) Wo + T^i is a better approximation than Lo + 

Example 5, Table I 

y"-\-ly'-\-y=0, OKx^co 



Solution required to match the accurate solution 






y{x) = Joix) - i No{x) 






at in&nity: 


■ 




1 +i -.-. 






\/irX 






1 + i ,v r . . / ^ .^ ■n\ 






"''-4V«'' r'^^ 'V''" 2). 


■ 




Tablk I 






X 


Jo - .Wo 


if'o 


if'o + Wl 


W!, + Wi + w. 


10 


- .2459 -/ .0557 


-.2468 -i .0526 


-.2460 -i .0557 






9 


- .0903 -1 


.2499 


-.0938 -i .2489 


-.0903 -/ .2500 






8 


.1717- 


.2235 


.1683 -i .2264 


.1717 -/" .2235 






7 


.3001 +1 


.0259 


.3009 +1 .0207 


.3001 +i .0260 






6 


.1506 +1 


.2882 


.1568 +i .2855 


.1507 +/ .2883 






5 


-.1776 +i 


.3085 


-.1704 +;■ .3135 


-.1777 +/ .3086 






4 


-.3975 +1 


.0169 


- .3979 +t .0291 


-.3973 +; .0169 






3 


-.2601 - 


.3769 


-.2765 -I .3684 


-.2601 -(: .3772 






2 


.2239 - 


.5104 


.1967 -/ .5288 


.2246 -i .5109 






1 


.7652 - 


.0883 


.7796 -i .1699 


.7683 -/ .0860 


.7651 -/ 


.0882 


0.8 


.8463 +1 


.0868 


.8920 -i .0130 


.8499 +/ .0916 


.8461 +i 


.0868 


0.6 


.9120 +( .3086 


1.0124 +i .1899 


.9152 +i .3183 


.9116 +/ 


.3084 
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For values of x less than 1 W2 was evaluated numerically. 
Example 6, Table II 

y{\) = 1, y'{\) = 0. 

The solution of this equation using Picard's method and the integraph 
has been described by Thornton C. Fry.^ We compare his results with those 
obtained by the wave perturbation metliod. The equation is first reduced 
to normal form by the substitution y — x~^'^ u, so that 



„" = -1 +^\u 



4a;2 



and we have /3 = ^s/S. Then 



,.1/2 



x"' Wo = cos d (x - 1) + -L sm 0{x - 1) 



2/3 



4/3 Ji \m2 



cos/3(» — 1) 



4- 2^ sin aiu 



-l)]si 



sin ^(.r — u) du. 



While W\ may be evaluated in terms of Ci and Si functions the values 
tabulated below were obtained by numerical integration. The values of 
the accurate solution 

y = 1.4034 /i(.i) - 0.3251 N^{x), 

and of the third and eighth Picard approximations, are copied from Fry's 
paper. 

Table II 



X 


y 


yi 


ya 


w„ 


iro + ifi 


1.0 


1.000 


1.000 


1.000 


1.000 


I.OOO 


1.2 


.998 


.998 


,998 


.990 


.998 


1.4 


.985 


.984 


.986 


.961 


.985 


1.6 


.956 


.951 


.955 


.913 


.956 


1.8 


.908 


.894 


.910 


.848 


.908 


2.0 


.842 


.809 


.844 


.769 


.842 


2.2 


.759 


.694 


.760 


.677 


.758 


2.4 


.659 


.548 


.661 


.575 


.659 


2.6 


.547 


.370 


.549 


.466 


.547 


2.8 


.425 


.156 


.427 


.352 


.427 


3.0 


.297 


-.096 


.300 


.236 


.300 
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